Origin of line tension for a Lennard-Jones nanodroplet 
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The existence and origin of line tension has remained controversial in literature. To address this issue we 
compute the shape of Lennard-Jones nanodrops using molecular dynamics and compare them to density 
functional theory in the approximation of the sharp kink interface. We show that the deviation from Young's 
law is very small and would correspond to a typical line tension length scale (defined as line tension divided by 
surface tension) similar to the molecular size and decreasing with Young's angle. We propose an alternative 
interpretation based on the geometry of the interface at the molecular scale. 



I. INTRODUCTION 

The development of microfluidics in the last decade has 
renewed the interest for a thermodynamical concept in- 
troduced by Gibbs in his pioneering article: line tension- . 
By analogy with surface tension, which is by definition 
the excess free energy per unit surface of an interface 
separating two phases, line tension is the excess free en- 
ergy per unit length of a contact line where three distinct 
phases coexist. The variation of a system free energy F 
therefore presents three contributions, a bulk contribu- 
tion when the volume V is varied, a surface contribution 
when any interface area 5^ is varied and a line contribu- 
tion when the contact line length L is varied: 

dF = PdV + lidSr + toIL . (1) 

i 

Here, we use a summation to indicate that one has to take 
all interfaces into account (liquid-solid, liquid-vapour, 
and solid- vapour). The stability of dcformable surfaces, 
such as a liquid-vapor or liquid-liquid interface, neces- 
sarily requires a positive surface tension. Although the 
shape of the contact line is deformable as well, the line 
tension cannot be inferred from a stability argument 2 -. 
In addition, there are conceptual problems defining line 
tension properly^. 

The simplest system in which a line tension effect may 
be observed is a liquid drop on a solid substrate, in partial 
wetting conditions. Consider a drop whose shape is a 
spherical cap characterised by its contact line radius R 
- seen from the top - and its contact angle 8. The drop 
volume is V — ^R 3 (2 — 3 cos 6 + cos 3 6*) , the liquid- 
vapour area Slv = 2nR 2 (l — cos6>), the solid- liquid area 
Ssl = 1" R 2 and the contact line length L = 2ttR. Here 
we defined R as the radius of curvature of the spherical 
cap: R = Rj sin 6*, see also Fig. [TJ When minimizing the 
free energy with respect to 9 at constant volume (PdV = 
0), one gets^: 

Isv-Isl T/7 t/7 

cos6>= — =cos6>y — , (2) 
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where 9y is the Young contact angle and 7sl, 7s v, 7 




FIG. 1. Schematic of two drops of same size, one cylindrical 
cap shaped (left) and the other spherical cap shaped (right). 
For small volumes, the contact angle 6 for the spherical cap is 
affected by line tension, while 9 is constant for the cylindrical 
cap. 



the solid-liquid, solid-vapour, and liquid-vapour surface 
tension, respectively. Note that ([2]) only holds for spher- 
ical cap-shaped droplets - the contact line of a cylinder- 
shaped drop has zero curvature, which means that the 
contact angle 9 is unaffected by line tension and indepen- 
dent of drop size. In this derivation, we did not take any 
interface curvature effects into account (such as Tolman- 
corrections on 7), if these become comparable in magni- 
tude to line tension the measured r from eq. [5] cannot 
be considered 'pure' line tension, but rather an appar- 
ent line tension 3 -^. From © one can see that when r is 
positive, drops will present a larger contact angle than 
Young's angle. 

Theoretical predictions on the strength of line tension 
are based on calculating the free energy (per unit length) 
associated with the contact line, using statistical mechan- 
ics (e.g. using density functional theory^) or a model 
based on interface displacement^. These analyses pre- 
dict the value of line tension to be in the range 10 -12 to 
10~ 10 J/m. Of particular interest is the behaviour near 
the wetting transition (6 — > 0), for which r can vanish or 
diverge depending on the details of the interaction-^— . 

A large amount of experimental work has been done to 
determine the magnitude of line tension. The most di- 
rect way is to measure the contact angle as a function of 
contact line curvature, and thus droplet sizei^— . Using 
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the modified Young-equation from ([2J r can then be cal- 
culated. Due to the small length scales involved for the 
measurement of r, the observed values for r vary greatly 
in magnitude: both negative and positive values as low as 
10~ n J/m and as high as 10 -5 J/m have been reported. 
The reason for the huge variation is that determining the 
contact angle is notoriously difficult due to contact an- 
gle hysteresis caused by surface inhomogeneities^. The 
slightest amount of surface inhomogeneities can cause a 
severe overcstimation of r. Indeed, contaminated sur- 
faces can even lead to an apparent change of sign of r 20 ' 21 . 
Historically, droplets were used to measure line tension. 
Recent developments on surface nanobubbles allowed to 
detect a similar size dependence of the contact-angle of 
the nanobubbles^—. Because of the difficulty of exact 
contact angle measurements at the required scale (1-100 
nm), alternative methods have been developed, for exam- 
ple by calculating the effective potential near the contact 
line by measuring the deviation of the liquid surface from 
a wedge shap o 26 i 27 . For reviews on experimental meth- 
ods and results see refs . 11 ' 16 . 

In this paper we adopt the usual experimental method 
to determine the line tension, by measuring 6 against 
1 /R, in a theoretical setting. We will perform this mea- 
surement for different equilibrium contact angles 9y, to 
study the dependence of the line tension length (£) on 
By. The tension length is defined as: 
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(3) 



Since we always find negative values of r for the Lennard- 
Jones drops studied in this paper, the minus sign is added 
to ensure that the tension length is always positive. 

We perform these measurements for both 3D (spheri- 
cal cap-shaped) and 2D (cylindrical-shaped) droplets, to 
compare similar sized droplets with and without contact 
line curvature. In the first part of this paper, we investi- 
gate line tension by means of molecular dynamics simu- 
lations of a Lennard- Jones droplet, which has the added 
advantage that no assumptions have to be made as is 
required for most analytical approaches. In this sense, 
these simulations are like an experiment, but with un- 
precedented accuracy and without surface inhomogene- 
ity. Since the expected tension length is of the order of 
the molecular size, the problem also has the right scale 
for molecular dynamics. Line tension has been observed 
in molecular dynamics studies befor o 28 ^ 29 , but a system- 
atic study has to our knowledge not been carried out . In 
the second part of the paper, we analyse the existence 
and origin of line tension using the density functional 
theory (DFT) in the approximation of the sharp kink 
interface. Finally, we will calculate the line tension us- 
ing a geometric interpretation based on missing bonds in 
a wedge-shaped interface. We show that the deviation 
from Young's law is very small and would correspond to 
a line tension of a fraction of the molecular size. 
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FIG. 2. Interaction potentials 4>(r) for the Lennard- Jones 
particles Q (solid line) and the potential p' 1 ' p' 2 ' g r (r)<f>(r) 
used for the DFT-calculations ([T7)l (dashed). Note that the 
DFT potential is regularized to account for vanishing g r {r) 
when r — ¥ 0. 



II. NANODROPS FROM MOLECULAR DYNAMICS 
A. Numerical setup 

Wc perform Molecular Dynamics (MD) simulations on 
nanodrops, using the Gromacs software package^. We 
simulate binary systems in which two types of particles 
exist: fluid particles that can move around either in the 
gas or liquid phase, and solid particles which are frozen on 
an fee lattice and constitute the solid substrate, Fig. (Ha) . 
The simulations are done in the NVT-ensemble, where 
the temperature is held at 300K using a thermostat, 
which is below the critical point for a Lennard- Jones fluid 
with the interaction strengths used. All particle interac- 
tions are defined by the Lennard- Jones (LJ) potential: 



(4) 



see also Fig. [2] Here, e,j is the interaction strength be- 
tween particles i and j and Oij the characteristic size of 
the molecules. This size is chosen the same for all inter- 
actions, oij — a. The potential function is truncated at a 
relatively large radius (r c = 5er) where (j> LJ is practically 
zero. The timestep is chosen at dt = a^J m/eu/200, with 
m the mass of the particles. The fluid particles are ini- 
tially positioned on an fcc-lattice near the substrate, but 
are free to move around and relax towards an equilibrium 
droplet shape (Fig. [3]). Periodic boundary conditions are 
present in all directions. To study the effect of line ten- 
sion we consider two different systems. In the "3D" case 
the dimensions of the system are chosen large enough to 
ensure that the droplet does not interact with itself, re- 
sulting in an isolated droplet with the shape of a spherical 
cap. In the "2D" case the system size in the x-direction 
(parallel to the substrate) is only 15a leading to an in- 
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finitely long cylindrical-cap shaped droplet. The small 
length is required to suppress the Rayleigh-instability, 
which is only effective at wavelengths A > 2wR. Now, 
as there is no contact line curvature, line tension has no 
effect. 

The wettability of the substrate (and thus, the equi- 
librium contact angle) is tuned through the parameter 
ratio thsl^LL-, where the subscript indicator S denotes 
the solid (fixed) particles and L the liquid particles. A 
higher value for cls/^ll results in a large attraction of 
fluid particles to the substrate, and thus a more wetting 
substrate. A range of contact angles can be explored in 
this way, as shown in table I. The obtained contact angles 
and compare well to those found in previous studie o 31 ' 32 . 
Depending on the size of the droplet, the effect of layering 
inside the liquid (Fig. [3]) limits the range where reliable 
contact angle measurements can be performed. In prac- 
tice, this limits the analysis to contact angles larger than 
approximately 70°. 

B. Cylindrical vs spherical caps 




FIG. 3. (a) Snapshots from Molecular Dynamics simula- 
tions of a 2D, cylindrical cap-shaped droplet (left), and a 3D, 
spherical cap-shaped droplet (right). The light spheres repre- 
sent the immobilized solid particles, forming the substrate to 
which the droplet attaches. The darker spheres represent the 
mobile fluid particles. The lines are a guide to the eye. Sev- 
eral periods of the 2D droplet are shown (periodic boundary 
conditions), causing the same particle to be printed multiple 
times. These drops were simulated using identical interaction 
parameters (|^ = | => 8y ~ 65°), and differ in shape only 
because of the difference in the periodic boundary conditions, 
(b) Isodensity contours measured using statistical averaging 
from the droplets shown in the top row. The contact angle 
and overall shape of the two drops is almost identical, re- 
quiring a precise measurement to observe the effect of line 
tension. 



Fig. Ela) shows the shape of two nano-drops (9y = 



65°) with similar radii R (as seen from the top), but 
with a different geometry. The cylindrical droplet on the 
left is formed in the quasi-2D system (several periods 
shown). The spherical cap shaped droplet on the right 
is simulated in a fully 3D system. Fig. [3jb) shows the 
isodensity contours from the same droplets. One observes 
that these cross-sectional shapes are already very similar, 
indicating that line tension is indeed a weak effect, even 
for such nanodrops. Wc will now perform careful and 
precise contact angle measurements in order to quantify 
line tension. 



C. Measurement of contact angle 

To perform precise contact angle measurements we first 
compute the density field by averaging over time and 
over space (translational or rotational symmetry). Dur- 
ing this averaging, we compensate for any center of mass 
motion of the droplet parallel to the substrate by mov- 
ing the droplet such that the center of mass is stationary 
throughout the averaging procedure. When the droplet 
has reached its equilibrium state (Fig. [3]), the density 
profiles are calculated by time-averaging over 1,000,000- 
10,000,000 time steps until the density field has con- 
verged. Using real world parameters for Argon as the 
fluid this would correspond to 2 - 20 nanoseconds. This 
leads to droplet shapes as shown in the bottom row of 
Fig. [3] The part of the droplet that is close to the sub- 
strate is subject to layering 3 ^: the density oscillates as a 
function of height. To avoid interference from this effect, 
we ignore this part of the droplet when determining the 
contact angle: we perform a circular fit through the top 
part of the spherical cap, and extrapolate towards the 
substrate (which is defined to be a/2 above the top row 
of substrate atoms) to find 9 and R. Fig.@Ja) shows these 
fits through some isodensity contours (p* — 0.3, 0.5, 0.7). 

This leads, however, to a new problem: which isoden- 
sity should one choose? As can readily be seen from 
Fig. SJa) it turns out that the width of the interface can- 
not be neglected, and choosing different isodensity con- 
tours results in different values for the 9 and R. To over- 
come this problem, we use the data from the cylindrical 
droplets to determine which isodensity contour to use. 
From a macroscopic perspective, the cylindrical caps are 
not affected by line tension (9 = 9y). It turns out that 
this property is obeyed by the Gibbs dividing surface at 
p* = 0.5, where p* is a parametrized version of the local 
density given by: 



Here, pl and pv are the bulk densities of the liquid and 
vapour phase, respectively. We note that although line 
tension does not affect cylindrical droplets, other curva- 
ture effects (such as the Tolman-correction on 7 and the 
effect of the increased Laplace-pressure on 7^;) do play 
a role. The baseline established by this methodology is 
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FIG. 4. Isodensity contours of a Lennard- Jones droplet, (a) A selection of isodensity contours (same drop as in Fig. Rq ps 9a, 
p* = 0.3, 0.5, 0.7) fitted by circles (dashed). The cirles turn out to be concentric which we use to collapse the contours into one 
single shape, (b) Density profile © fitted to measured p*{r), where we can see that the interface is several molecular diameters 
thick. With this fit we can shift all contours by using (B). (c) Rescaled iso-contours nicely collapse on a single curve, which 
allows us to define the interface in a precise way. 



therefore not based on 'pure' line tension, but rather an 
apparent line tension in which all these effects are com- 
bined. This is in line with previous experimental work, 
where this distinction could also not be made. 

To improve statistics, we have also used the remaining 
isodensity contours. This can be done since the den- 
sity profile across the interface is accurately fitted by 
(Fig. lib)): 

"* = K 1 + te " h (^))' TO 

where Rq is the point where the density is halfway be- 
tween the liquid value and the vapour value (p* = 0.5). 
w is a fit parameter that defines the width of the liquid- 
vapour transition. Since the circular fits are concentric 
(they share a common center point, C in Fig. Ha)) we 
can easily transform any isodensity contour to the refer- 
ence contour. For this, we calculate the radial distance 
from the contour towards the reference contour using (|6|) . 
The result of this transformation for the spherical droplet 
from Fig. [3] is shown in Fig. He), where we see that the 
contour shapes indeed collapse and can now all be used 
to determine the contact angle. The spread of these val- 
ues for different p* are used to determine the error of the 
measurements. 
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FIG. 5. cos# vs aRr 1 for cylindrical (open symbols) and 
spherical (filled symbols) drops for cos8y = —0.60 (« 127°). 
The dashed lines are linear fits through the data points, and 
the solid lines are the solutions obtained using DFT described 
in section IITT1 The top two (red) lines represent the 3D-data, 
whereas the bottom two (green) lines the 2D-data. The di- 
amond at 1/R — indicates Young's law, calculated inde- 
pendently by determining the surface tensions of the three 
interfaces: 7, fsL, ysv- The difference between the slopes of 
the 2D and 3D fits quantifies the tension length i. Note that 
for this particular equilibrium contact angle, MD and DFT 
agree quantitatively on the tension length. 



D. Results: tension length 

Figure [5] shows the relation between the contact an- 
gle defined as described above and the drop radius for 
different sized droplets. Young's angle By was indepen- 
dently calculated from independent measurements of the 
surface tensions^ of planar interfaces, under the same 
simulation conditions as the droplet simulations. This 



is shown as the diamond symbol at 1/R, = 0. One can 
observe that the contact angle does not present any vari- 
ation with drop size in the cylindrical cap case, which 
is consistent with the macroscopic picture. By contrast, 
the spherical drops exhibit a decreasing contact angle for 
small radii (large i? _1 ), which according to @ is con- 
sistent with a negative line tension r. A negative value 
of r means that the contact line has the tendency to ex- 



5 




FIG. 6. Tension length I vs Oy for spherical drops. Square 
symbols are the molecular dynamics results, triangle sym- 
bols are the results from the self-consistent Density Functional 
Theory model discussed in section IIII Al These data points 
were acquired by measuring the contact angle for different 
drop sizes, meaning they represent an 'apparent' line tension. 
The solid and dashed lines also result from DFT, assuming 
a wedge-shaped geometry near the contact line: eqs. (|25|l . 
(I26|l . (|27p . For the self-consistent DFT data the characteristic 
lengths are determined analytically: C,ll = Cls ~ no /A. The 
resulting curve is the solid line. The characteristic lengths for 
the MD data are acquired by fitting: (ll = 3.5 a, (ls = 0, 
represented by the dashed line. 
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TABLE I. MD results as shown in Fig. [5] The ratio e LS /e LL 
was varied to obtain the tension length I for different equilib- 
rium contact angles By. #dft is the contact angle resulting 
from the interaction ratio according to the DFT-model de- 
scribed in Sect. IIHI 

pand - a larger contact line length leads to a decrease in 
9 under the constant volume constraint. The solid line 
corresponds to the density functional theory in the sharp 
kink approximation that will be discussed below. 

The difference between the slopes of the 2D- and 3D- 
fits in Fig. [3] (dashed lines) is equal to the tension length, 
£ = — r/7, which is defined to be positive for negative 
values of r (see ©). For this equilibrium contact angle 
(&Y = 127°) we find £ = 0.36cr. Now, by varying the 
interaction ratio cls/^ll we measure £ for varying 9y. 
The result is shown in tableland in Fig.|6]by the square 
symbols. Whatever 9y, the tension length turns out to 
be positive (so r is always negative) and very small - of 
the order of the atomic size a. The tension length recov- 



ered from the MD simulations is a decreasing function of 
9y , indicating that the effect is stronger when the wedge 
formed by the liquid in the vicinity of the contact line is 
sharp. The other curves in Fig. [5] result from DFT, and 
will be discussed in the following sections. 



III. ORIGIN OF LINE TENSION EFFECT 

In this section we study line tension in the framework 
of Density Functional Theory using the sharp kink ap- 
proximation. Once more, the strategy is to determine the 
equilibrium shapes of 2D and 3D drops and to compare 
their contact angles. Starting from the basic equations of 
DFT we first motivate the form of the free energy func- 
tional in Sec. IIII Al Some of the assumptions are directly 
tested using Molecular Dynamics simulations. We then 
derive the equilibrium condition for the capillary pres- 
sure (Sec. IIII Bp and describe the numerical scheme that 
was used to solve the equilibrium shapes of the drops 
fSec. IIIICl) . The numerical results are presented and in- 
terpreted in detail in Sees. IIII Dl and IIII El 



A. Density Functional Theory in the sharp kink 
approximation 

The primary idea of Density Functional Theory (DFT) 
is to express the grand potential ft = U — TS — [iN = 
F - fiN BjS ci functional of the particle density p and 
to perform a functional minimisation for a given \x and 
T. For an ideal gas, the free energy functional is known 
explicitly 

F id [p] = kT jp [ln(M 3 ) - 1] dr , (7) 

but this is not the case for general liquids. Let us denote 
0(^1 j ^2) or (jj(r) as the additive pair potential between 
particles at r\ and r 2 with distance r = \f 2 — ri|. From 
a g rand canonical averaging, one can show (see e.g. refs. 
HJandii) 

7T7^ - \ = xP^Hn,^) = lp(fi)p(f 2 )g(fi,f 2 ) , (8) 
d0(ri,r 2 ) 2 2 

where is the two-body density distribution function 
and g the pair correlation function. This relation can be 
used to construct the free energy for non-ideal systems. 
Introducing a coupling parameter A in front of the inter- 
action, the free energy can be constructed by integration 
as: 

F[p] = F id [p] + 

^ J d\ j d?i j df 2 p(r 1 )p(f 2 )gx(ri,r 2 )(l)(\r f 2-n\) ■ 

(9) 
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of the solid and liquid, we find 



FIG. 7. Schematic representation of the integration variables 
and domains from (|I1[) . Both the liquid-solid (left) and the 
liquid-liquid (right) interactions are integrated over their re- 
spective volumes (the liquid cap C and the solid substrate <S) 
to obtain the total free energy F . 



Here g\ is the pair correlation function in a system of 
same geometry and same volume, for which the interac- 
tion is X<p(r). 

Although exact, this expression cannot be used as it 
is, as g\ is not known. For a practical approximation of 
the energy functional, one can separate the thermody- 
namic non-ideality in contributions due to attractive and 
repulsive components of the intermolecular potential. As 
the repulsive forces have a very short range, their effect 
is mainly local. Using the Local Density Approximation, 
the repulsive contribution can be estimated from the Hcl- 
moltz energy density f r (p) in a uniform system of den- 
sity p at temperature T, composed of purely repulsive 
molecules. The attractive Van der Waals interactions, 
4> a tt, can then be treated as a perturbation, assuming 
that the pair correlation function remains mostly that of 
the purely repulsive reference system, g r (r*i, r^)- The free 
energy then reads: 



F[ P ] = / f r (p)df + 



- I dri I dr 2 p{fx)p{f 2 )g r {ri,r 2 )(j) at t{\r 2 - ri\) 



(10) 



F = Mpl) j^dr + f r (p s ) J dr 

+ zfif dft J dr 2 g r (\r 2 - ri\)<j> LL (\r 2 -f\\) 

+ PlPs J dri J dr 2 g r {\r 2 - n |)0Ls(|r 2 - fi|) , 

(11) 

where C and S are the liquid and solid domains respec- 
tively, see also Fig. [7] <j)LL and (/)ls denote the attractive 
parts of the respective interactions. The drop shapes are 
found by minimizing this free energy with respect to the 
shape of the liquid domain C. 

Before proceeding, it is instructive to discuss the ap- 
proximations underlying pT|) in the light of the Molecu- 
lar Dynamics simulations of the Lennard- Jones droplets. 
First, the assumption that the pair-correlation function is 
homogeneous in space ignores the layering near the solid 
wall (cf. Fig. [3]). This can induce significant corrections 
to the estimated free energy. Second, the Local Density 
Approximation of the short-range repulsive forces gives 
rise to isotropic repulsive interactions, while the attrac- 
tive interactions will become anisotropic in the vicinity 
of an interface. If this is indeed the case, the surface 
tension (and line tension) effects mostly result from the 
attractive component of the interaction. We test the va- 
lidity of this hypothesis in the Molecular Dynamics sim- 
ulations by measuring the anisotropy of the stress tensor 
in the vicinity of the liquid- vapour interface. We define a 
cumulative stress-tensor a aa (R*) that incorporates only 
the interactions with a bond length smaller than R* : 



'{R*) 
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(12) 



The true stress in the system is recovered when R* = oo, 
for which all interactions are taken into account. Here, 
rrii, Vi are the mass and velocity of particle i, respec- 
tively, and fij and r,j are the force and displacement 
vector between particles i and j. With this, we quan- 
tify the anisotropy from the difference between the stress 
components tangential (T) and normal (./V) to the inter- 
face, as 



To end up with a numerically tractable scheme, we 
make a final approximation that the density profile across 
the interface is mostly independent of the geometry. 
Defining the position of the interface e.g. by the iso- 
density p* — 1/2, the integrals in (flU)) can be approxi- 
mated by assuming that the density is uniform in both 
phases^.This so-called 'sharp kink approximation' ne- 
glects the thickness of the diffuse interface. Thermal ef- 
fects are implicitly taken into account, since f r , g r and 
the liquid density depend on temperature. In this ap- 
proximation, the free energy becomes an explicit func- 
tional of the shape of liquid, solid, and vapor domains. 
Since the vapor density is neglible with respect to that 
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(13) 



Figure [5] shows the anisotropy A as a function of R*. The 
dashed line indicates the transition from the repulsive 
(r < 2 1 / 6 cr) to the attractive domain (r > 2 1 / f V). The 
figure clearly shows that the majority of the anisotropy in 
the liquid-vapour interface is due to the attractive inter- 
action, while the repulsive interaction accounts for about 
20% of the anisotropy. This indeed justifies a local den- 
sity approximation for the repulsion, although one can 
expect quantitative differences with Molecular Dynam- 
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FIG. 8. Stress anisotropy yl for bond lengths smaller than 
R* , see (|13p . The measurement was done in a slab of height 
a/3 within the liquid- vapour interface in a molecular dynam- 
ics simulation. The dashed line indicates the minimum of 
the Lennard- Jones potential at R* = 2 1 ^ 6 cr and marks the 
separation between the attractive and repulsive bonds. The 
majority of the anisotropy comes from attraction. 



B. Capillary pressure 

The equilibrium shape of liquid drops can be obtained 
by minimizing the free energy F at constant volume V. 
This can be done by variation of (| 1 1 1) with respect to the 
drop shape C under the constraint of constant volume. 
The resulting equilibrium condition is a constant poten- 
tial energy density IT along the free surfac o 6 ' 36 ' 37 . This 
potential can be interpreted as the capillary pressure and 
can be decomposed into a liquid-liquid and a solid-liquid 
contribution, as II = Ull + ^ls- The former can be 
written as 

Il LL (r) = -H° LL + PL J df ft-fl? - A)4>ll{\? - r\) , 

(14) 

where we subtracted n^ L , the interaction due to a semi- 
infinite volume of liquid. The solid-liquid contribution 
follows from the interaction due to the semi-infinite vol- 
ume of solid 

IT-ls (r) = PS [ d?g r (\7 - r\)<b LS (\7 - r\) . (15) 
Js 

The equilibrium condition is thus that 

n (**) = U LL (r*) + n LS (r*) = constant, (16) 

where r* denotes an arbitrary position at the liquid- vapor 
interface. 

Note that the capillary pressure IT depends on the 
shape of the liquid, through the domain of integrations, 
and thus reflects the effect of the interface geometry 
on the free energy. It implicitly contains the Laplace 
pressure, which is the capillary pressure associated to a 
macroscopic curvature of the interface, and the disjoining 



pressure, which is the capillary pressure in the case of a 
microscopic film. 

C. Shape relaxation 

We compute droplet shapes for a pair interaction con- 
sistent with the long-range van der Waals interaction 
used in the Molecular Dynamics simulations (see Fig. : 

PiP 3 9r{r)(t)ij{r) = lJ 3 , (17) 

[o + r z ) 

where represents the strength of the interaction be- 
tween molecule i and j. Comparing to the van der Waals 
interaction of ((4]), one finds Cjj = 4piPja 6 eij. For mathe- 
matical convenience we have chosen a simple rcgulariza- 
tion around r = a, which represents the effective size of 
the short-range repulsion. Let us note that the poten- 
tial of equation (17) does not lead to the formation of 
a precursor film. Namely, the corresponding energy per 
unit surface for a flat film is a monotonic function of the 
thickness h, with a prefactor depending on the spread- 
ing parameter. For partial wetting, the system tends to 
zero thickness rather than to a precursor film of finite 
thickness. We have tested that other similar choices like 
g(r) = for r < a and g(r) = 1 for r > a leads to 
quantitatively similar results (see ref. I37T ). The surface 
tensions corresponding to (|17[) can be computed as 

7 = TTc u /8a 2 , (18) 

and 

7 + Isv - J si = Trcis/Aa 2 . (19) 

By choosing identical functional forms for both interac- 
tions, one simply has 3 ^ 

COS0V- = COS 9d ft = 2— - 1. (20) 

Similar to the Molecular Dynamics simulations, we 
compute the equilibrium shapes of nanodrops in both the 
2D configuration (cylindrical caps) and 3D configuration 
(spherical caps). The drop shapes are parameterized by 
r(a) as shown in Fig. [5] - polar coordinates are used 
to allow contact angles larger than 7r/2. We numerically 
determine the equilibrium shape of the drop by an it- 
erative algorithm that tends to a constant IT(a) along 
the interface. The initial shape is taken as a spherical 
cap with 0y according to (|20l) . This is shown in Fig. [9] 
by dashed lines. The corresponding potential 11(a) is 
uniform except within a few molecular scales from the 
contact line, where the influence of the solid plays a role. 
We iteratively construct drop shapes r*(a) according to 
r t+1 (a) = r*(a) + X t (n t (a) — (7r*) Q ), while keeping the 
volume constant. Here, 7r*(a) is the capillary pressure at 
angle a during iteration t. (tx } is the space-averaged 



FIG. 9. Typical result of the DFT model in the sharp kink 
approximation, (a) Surface potentials of an axisymmetric 
drop with 6y = 127°, a measured angle 9 = 116° and ra- 
dius R fa 2cr. a = at the contact line. The dashed line 
(top) represents ft;;, the bottom line represents II S ;, and the 
lines in the middle indicate the sum of the two. Here, the 
dashed and solid lines indicate the potential energy density 
of the droplet in its initial shape, and its equilibrium (final) 
shape, respectively, (b) Initial and final drop profiles shown 
by the dashed and solid lines, respectively. The initial profile 
is a spherical cap shaped drop with 8 = 9y = 127° . 
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FIG. 11. Cosine of the equilibrium contact angle against 1/R, 
for a 3D drop By = 28° as shown in the insets for two different 
sizes R. At sufficiently small radius we observe a saturation of 
cos(#) = 1, approaching a perfectly wetting drop (top inset). 
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FIG. 10. Cosine of the equilibrium contact angle against 1/R, 
for 2D drops (triangles) and 3D drops (squares). The corre- 
sponding Young's angle is 9y — 65° . The slope of the curve 
near 1/R = can be attributed to line tension for the 3D 
drops, while it is zero for 2D drops. Note that both curves 
exhibit a significant 1 / R 2 contribution for smaller drop sizes. 
This contribution was not recovered from the molecular dy- 
namics simulations, since the radius of the droplets was not 
small enough: R > la. Smaller droplets would not allow for 
spherical cap fitting because the droplet size becomes similar 
to the particle size. 



potential at the interface during iteration t. The param- 
eter A* is selected such that the variance of the potential 
is minimised at each step. After a few hundred steps, the 
shape converges and yields 11(a) that indeed is constant 
within numerical precision. Note that for the potential 
studied here, no precursor film is formed. 

The shape r(a) and the capillary pressure 11(a) of a 
small drop are plotted as solid lines in Fig. |9] Away from 
the contact line the drop is a spherical cap, but a sig- 
nificant deviation can be observed near the contact line. 



The drop has spread with respect to the initial shape, 
resulting into a lower contact angle than 6y ■ Once more, 
this is consistent with a negative value of the line tension 
r. Far from the contact line, the capillary pressure is 
dominated by the II;; term. The corresponding value is 
simply the expected Laplace pressure 2^/R, where R is 
the radius of curvature the drop. 



D. Results 

In Fig. [TU] we compare the contact angles of 3D drops 
(squares) and 2D drops (triangles), as a function of in- 
verse drop radius 1/R. In both cases the interactions 
were identical, corresponding to Oy = 65°. For large 
drops there is indeed a difference that can be attributed 
to line tension: the slope at 1/R — > is finite for 3D 
drops while it vanishes in the 2D case. Interestingly, how- 
ever, there remains a 1/R 2 contribution for both types 
of drops. The two data sets are accurately fitted by 
parabola, with equal prefactors for the quadratic term. 
This suggests that the effect of line tension in (J5|) can be 
seen as the leading order contribution of an expansion 
in a/R, and is only valid for relatively large drops. In 
particular, ^ must break down when cos 8 « 1. This 
is illustrated by Fig. [TT] showing a saturation of the con- 
tact angle to 8 rj for very small drops. This effect is 
of course most pronounced for drops that already have a 
small Young's angle 9y . For such small drops, the range 
over which one observes a 1/R behaviour is very small 
and the main size effect is to induce a wetting transition. 

We are now in the position to make a comparison of 
the DFT model with the molecular dynamics simulations 
presented in Sec. HU The solid line in Fig. [5] represents 
the contact angles for 3D drops of varying sizes as ob- 
tained from the numerical DFT calculation. We took the 
same Young's angle as obtained in the molecular dynam- 
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ics simulations, i.e. By = 127°. The trends of DFT and 
molecular dynamics are very similar, clearly showing a 
decrease in contact angle for decreasing drop radius. For 
both cases line tension is thus negative and has a similar 
magnitude in units of a. 

Finally, we determined the tension length t from the 
slope near 1/R — > 0, for a broad range of By. The re- 
sults are reported as triangles in Fig. [51 The value of £ 
vanishes both for and 180° and presents a maximum 
around 90°. Beside the sign and the order of magni- 
tude, the behaviour is thus qualitatively different from 
the molecular dynamics results. This difference is most 
pronounced at small By (Fig. [6]). 



E. Geometric interpretation of line tension 

Within our DFT model, the dependence of i on contact 
angle By can be accurately described from a geometric 
argument (solid line in Fig. [5]). We separate the free en- 
ergy in volumic, surfacic and linear contributions, as 
F = PV + jiSi + tL. By assuming the liquid domain 
to be wedge-shaped, it is indeed possible to explicitly sep- 
arate the domains of integration in in bulk, surface, 
and line contributions: 

F = PV + J2 liSi 

i 

+ t;Pl dn dr 2 g r {\f2 - fi\)(j)LL{\r2 - n\) 
1 Jc[ Jc> 2 



+ PlPs I dfi / dr 2 g r {\r2 - ri\)4> L s{\r 2 - r\\) 
J c Js' 



(21) 



Here the integration domains £[, C 2 , and S' are those 
represented in Fig. Upland the Appendix. Note that such 
a decomposition is uniquely defined in the sharp-kink 
approximation, while this is no longer the case for inho- 
mogeneous density profiles. 

From (f2"Tj) one sees directly that line tension has two 
contributions, due to liquid-liquid interactions tll , and 
due to liquid-solid interactions tls- These can be com- 
puted as follows. The integration domains C' l7 C 2 , 
and S' arc bordered by straight lines passing through 
the contact line so that they do not present a character- 
istic scale. Therefore, both tll and tls can be written as 
products of a characteristic length (that does not depend 
on 6), and a function of 9 (that does not depend on the 
potentials </)ll and 4>ls)- It turns out that the lengths 
can be expressed in terms of the liquid-liquid and solid- 
liquid disjoining pressures U^(h) and IL^g(h). The dis- 
joining pressure is the energy per unit liquid volume at a 
distance h from a flat semi-infinite zone of liquid or solid 
(see the integration domain S' in Fig. [T2"]) . The surface 
tensions, already computed in (|18I19|) . can be expressed 




FIG. 12. Integration domains of the free energy (|11[) that 
contribute to line tension due to liquid-liquid interactions tll 
(left) and solid-liquid interactions tls (right). Assuming the 
contact line region to be perfectly wedge-shaped the total free 
energy F minus the volumic and surfacic contributions results 
in a residual energy, which can be attributed to line tension. 
We show in the Appendix how the line tension contribution 
can be isolated from the liquid-liquid and solid-liquid inter- 
actions. Then, by calculating the free energy associated with 
these integration domains (Eq. ([21])) one directly finds the 
line tension due to liquid-liquid interactions and solid-liquid 
interactions: eqs. (f25j, ([26]), and ([27)l . 



as the integrals of these quantities: 



Ilf2(h)dh = 2 1 (22) 



nff (h) dh = j + lsv - 7si = 7 (1 + cos By) (23) 



The characteristic lengths Cll and Cls that appear in 
the calculation of the line tension turn out to be the first 
moment of the disjoining pressure: 



Cll 



and Cls 



J zU Ls( z )dz 
JILf§(z)dz 



(24) 

Following the interpretation of surface tension as a force 
per unit length, Qll and Qls are the "moment arms" of 
these forces. Within our DFT model, the liquid- liquid 
and solid-liquid potentials have the same shape so that 
these two lengths are equal, (ls = Cll = 7rcr/4. 
The line tension tls follows as: 



tls — ClsI 



(1 + COs6»y) 

tan 8 



(25) 



This contribution is positive for < B < ir/2 and changes 
sign at 8 = ir/2. The prcfactor (1 + cos 6>y) is not of geo- 
metric origin, but stems from the strength of the liquid- 
solid interaction c^s- A similar result for t l $ was pre- 
viously obtained in ref. HU, but this work omitted the 
contribution due to liquid-liquid interactions, tll, which 
is crucial to describe our numerical DFT results. The 
contribution due to liquid-liquid interactions is negative 
for all angles. In the limit of small angles, t l l diverges 
as 



tll 



-Cll 1 



tan6> 



(26) 
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As the angle 9 tends to tt, tll vanishes as 

tll = -^-Clli(tt-0) 2 (27) 

In between, we have determined the ratio tll/Qll by 
numerical integration. 

Adding the two contributions tll and tls, wc obtain 
the solid line in Fig. HI which indeed closely follows the 
full numerical simulations obtained from the spherical 
cap measurements. Note that both tll and tls scale 
as 1/9 for small angles, but the diverging contributions 
balance exactly. This is a consequence of having identical 
values for the moment arms, i.e. = Cls, resulting in 
a vanishing line tension for small 9. Of course, this will 
not be the case in general, where wc expect one of the 
contributions to dominate. 



IV. DISCUSSION 

We theoretically investigated the effect of line tension 
by studying the contact angles of Lcnnard- Jones droplets 
of varying sizes. The equilibrium shapes of nanodrops 
were determined using two methods: Molecular Dynam- 
ics (MD) and Density Functional Theory (DFT). For 3D 
drops wc found a size-dependent contact angle consistent 
with while the contact angle was nearly constant for 
2D drops. DFT in the employed approximation does not 
fully reproduce the MD simulations, but it does capture 
the main physics. In particular, DFT gives the correct 
(negative) sign and order of magnitude of t, and also 
captures the dependence on wettability for large contact 
angles. Obiously, the exact numerical values resulting 
from the DFT calculation depend on our specific choice 
of the potential (|T7|) . Note however, that both the re- 
covered trend and the orders of magnitude for £ are a 
general result, independent of the specific choice of the 
potential. The only exception is the limit of the wetting 
transition, 0—5-0, which is known to depend on details 
of the interactional 2 *^. 

In addition, we identified a simple geometric interpre- 
tation of line tension. Molecules inside a liquid wedge 
interact with a larger number of surrounding molecules 
than estimated from surface tension, which is based on 
an infinite half space of liquid. Hence the negative sign 
of line tension. The wedge shape of the liquid is indeed 
a good approximation of the liquid geometry for large 
contact angles and yields a very accurate prediction for 
r in the DFT case. This is remarkable, since these DFT- 
measurcmcnts did not discriminate between line tension 
and other curvature effects, suggesting that line tension 
is the dominant mechanism for the size dependence of 
the contact angle. Once more, the behavior for small 
contact angles is sensitive to details of the interaction: 
it depends on the "moment arm" of the surface tensions, 
characterized by the length scales Qll and £ls- We spec- 
ulate that the layering effect near the substrate in MD 



substantially reduces the moment arm (ls for the liquid- 
solid interaction. This would explain the discrepancy 
with DFT. Indeed, the MD data can be described by 
the wedge-approximation of t^l by fitting the moment 
arms to C,ll = 3.5er and Cls = 0. It would be interesting 
to further investigate this matter. 

Although we were able to observe the variation of con- 
tact angle with drop size, the effect is only noticeable 
for very small, nano-scalc drops. Taking a = 0.34 nm 
and 7 = 0.017 J/m 2 , our results correspond to line ten- 
sion in the range r = 10~ 12 — 10~ n J/m (depending 
on the wettability). This is consistent with theoretical 
predictions as well as with recent experiments 2 ^. Note, 
however, that much larger experimental values for r have 
also been reported^—. Resolving this issue is particu- 
larly important for surface nanobubble a 22 i 25 ' 39 ' 40 , typi- 
cally lOOnm wide, whose stability was suggested to rely 
on an effective line tension^. 
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Appendix A: Geometric interpretion of line tension 

Figure [T3] shows by illustration how the free energy 
associated with the liquid-liquid interactions of a wedge- 
shaped liquid interface sitting on a solid can be decom- 
posed in its bulk, surface and contact line contributions. 
Fll shows the total free energy of the liquid-liquid inter- 
actions, which is the interaction of the liquid in the wedge 
with itself. For clarity, we separated the two domains in 
the first row spatially, but in reality they of course over- 
lap since they are the same volume of liquid. First, we 
decompose the integral domain in the bulk energy con- 
tribution (wedge shape ® infinite volume): -Ftmik- The 
surplus that has to be subtracted is shown in the top right 
of Fig. [T3l because one has to compensate for the areas 
where no liquid is present. From this surplus we extract 
the surface contributions. Note that the liquid wedge has 
two surfaces: the liquid- vapour interface and the liquid- 
solid interface, which are both represented by integration 
of the wedge (dotted area) with an infinite half-space, re- 
sulting in the total surface energy term. The third row 
shows what remains and is by definition (Eq. JT])) the line 
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FIG. 13. Integration domain for the liquid-liquid interaction energy decomposed in the bulk, surface, and line components. 
The dotted and striped regions represent the domains of integration for the variables df\ and dr*2, respectively, in the liquid- 
liquid term of (|21|l . The remainder after subtracting the bulk energy (fbuik) and surface energy (Fsurface) is the free energy 
associated with liquid- liquid line tension (tll)- Note that the two line-tension contributions shown here can be combined into 
the integration domains shown in Fig. [T2j 



tension. These integration domains can be simplified and 
merged into the one shown in Fig. [12] (left). 

To compute tls we follow a similar route. Fig. [H] 
(left) shows the integration domain for Fls for a liquid 
wedge (dotted) in contact with a solid (striped). The 
right two panels directly give the decomposition into the 
surfacic component (solid half-space ® liquid half-space, 
over the solid-liquid interface), and the remainder which 
is the line tension component (tls)- There is no bulk 
energy term since we are dealing with two separate (and 
spatially separated) phases. 
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FIG. 14. Integration domain for the liquid-solid interaction energy (left). This integration domain can be decomposed in the 
corresponding surface energy contribution (-Fsurfacc) and the free energy associated with liquid-solid line tension (tls)- The 
dotted region represents the integration variable dfi in the liquid-solid term of (|21|l and the double striped region the integration 
variable dr-z in the same equation. 



